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Summary 

During  the  tenure  of  this  research  grant  the  principal  investigator  (PI)  worked  in  the  areas  of 
random  iterated  function  systems  (IFS)  and  computer  image  generation.  Random  IFS  afford  a 
powerful  new  technique  for  generating  images.  The  image  is  constructed  as  the  attractor  for  a  random 
dynamical  system  in  R3  or  R3,  typically  a  Markov  process  or  some  spin-off  of  such  a  process.  The 
image  generation  algorithm  proceeds  by  simply  generating  the  appropriate  random  dynamical  system 
(via  a  random  number  generator)  and  plotting  the  points  along  a  single  orbit.  The  random  dynamics 
evolve  through  compositions  of  independent  identically  distributed  (i.i.d.)  affine  transformations. 

The  images  obtained  in  this  way  include  fractals,  landscape,  natural  scenes,  simple  geometrical 
shapes,  smooth  curves  and  surfaces,  rough  curves  and  surfaces,  wavelets,  B-splines,  Bezier  curves  and 
interpolants.  In  fact  this  technique  produces  one  of  the  most  efficient  and  parallelizable  algorithms  for 
high-speed  curve  and  surface  generation. 

Key  questions  are  how  to  find  the  appropriate  random  dynamics  to  generate  a  given  target 
image  (the  “encoding  problem”),  and  how  to  adapt  the  algorithm  so  as  to  produce  different  classes  of 
images.  There  were  several  research  accomplishments  in  these  directions  under  this  grant.  They  are 
listed  here,  and  elaborated  on  below  in  the  sections  which  follow. 

I)  Encoding  via  the  Collage  Property: 

Together  with  J.-P.  Vidal  (Ph.D.  student,  Carnegie  Mellon  University)  the  PI  developed  an 
automated  encoding  scheme.  This  encoder  works  by  sequentially  coilaging  a  given  polygonalized  target 
image  by  affine  copies  of  itself.  It  essentially  solves  a  puzzle  by  moving  and  adjusting  the  pieces  locally 
in  an  iterative  fashion,  so  that  they  efficiently  collage  the  target.  (This  is  not  quite  the  same  as  a 
puzzle  —  since  the  pieces  are  allowed  to  overlap  here.) 

II)  Encoding  via  Convexity: 

Together  with  J.-P.  Vidal  the  PI  discovered  an  important  property  satisfied  by  the  extreme 
points  of  (the  convex  hull  of)  the  attractor  for  an  IFS.  Namely,  every  extreme  point  is  the  image  of 
some  extreme  point  under  one  of  the  affine  maps.  Equivalently,  the  set  of  extreme  points  is  invariant 
under  the  inverse  dynamics.  Thus  by  examining  the  extreme  points  of  the  attractor,  one  can  extract 
information  about  the  affine  maps  which  generate  it.  This  extreme  point  phenomena  was  observed  by 
the  PI  to  hold  as  well  for  more  general  iterated  function  systems  than  affine  —  for  example  Julia  sets 
for  iterates  of  complex  polynomials. 
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III)  Encoding  via  Discrete  Optimization: 

Together  with  M.  Perrugia  (Ph.D.  student,  Carnegie  Mellon  University)  the  PI  studied  the 
discrete  Markov  chain  obtained  by  truncating  the  continuous  process  in  R3  to  the  centers  of  the  pixels 
in  which  it  lands.  It  was  shown  that  as  pixel  dimension  decreases  to  zero,  stationary  distributions  for 
the  discrete  truncated  chains  converge  weakly  to  that  of  the  continuous  IFS  process.  This  then  enables 
one  to  approach  the  encoding  problem  through  discrete  optimization. 

IV)  Encoding  via  the  Radon  Transform: 

The  Radon  transform  of  the  invariant  measure  for  an  IFS  satisfies  a  difference  equation,  with 
difference  parameters  related  to  the  affine  maps  of  the  IFS.  Knowing  the  target  image  as  a  color  or 
grey-level  image  amounts  to  knowing  the  invariant  measure,  and  from  this  one  can  obtain  its  Radon 
transform  via  plane  wave  scans.  So  for  this  approach  encoding  becomes  the  inverse  problem  of 
recovering  the  parameters  of  a  difference  equation  from  measurements  of  its  solution  (the  Radon 
transform). 

V)  Mixing  Images: 

Together  with  M.  Barnsley  (Iterated  Systems,  Inc.)  and  M.  Soner  (Carnegie  Mellon 
University),  the  PI  developed  a  framework  for  mixing  image  textures.  For  example  one  could  start 
with  a  given  leaf  texture  and  tree  texture,  and  then  create  a  tree  with  those  leaves  growing  on  its 
branches.  The  mixing  involves  constructing  a  stochastic  process  which  jumps  between  various  Markov 
processes,  transferring  states  from  one  to  the  other  via  a  graph  structure.  The  nodes  of  this  graph 
represent  simple  (un-raixed)  IFS  processes,  and  the  edges  represent  crossovers,  whereby  a  point  in  the 
orbit  of  one  process  gets  replaced  by  the  corresponding  point  in  the  orbit  of  an  acfjacent  process. 

VI)  Real-Time  Animation  and  Vectorization: 

IFS  images  can  be  animated  by  continuously  varying  the  IFS  parameters  (the  coefficients  of 
the  affine  maps  and  the  probabilities  assigned  to  them).  The  IFS  algorithm  is  highly  parallel,  and 
animation  is  also  parallel  since  there  are  no  dependencies  between  one  image  and  the  next.  (That  is, 
they  can  all  be  generated  simultaneously.)  The  PI  used  the  vectorization  and  multi-processor 
capabilities  of  the  CRAY  Y-MP/832  at  the  Pittsburgh  Supercomputing  Center  to  optimize  IFS 
animation,  thereby  generating  the  frames  in  real-time  (30  frames/sec.).  This  work  also  demonstrates 
the  continuous  dependence  of  the  image  on  the  IFS  parameters,  as  continuous  variation  of  parameters 
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leads  to  a  continuous  flow  of  images. 

VII)  Wavelets  and  Curve  Generation 

Most  recently  the  PI  discovered  how  to  set  up  a  random  IFS  to  generate  wavelets  and  various 
curves  arising  from  sub-division  methods.  These  include  B-splines,  Bezier  curves  and  line  averaging 
interpolants.  In  fact  for  a  whole  variety  of  computer  image  generation  algorithms  involving  recursive 
tree  traversal,  the  PI  has  shown  how  one  can  construct  a  simple  random  process  which  generates  the 
same  image.  The  image  then  evolves  as  the  orbit  of  a  single  trajectory.  This  involves  products  of  i.i.d. 
row  stochastic  matrices,  as  described  below. 


Introduction 


The  basic  IFS  image  generation  algorithm  is  illustrated  in  Figure  1.  The  leaf  is  generated  as 

follows.  Pick  any  point  X0  €  R2.  There  are  four  affine  transformations  T:  x  Ai  +  b  listed  on  the 

top  of  this  Fig.,  and  four  probabilities  p<  underneath  them.  Choose  one  of  the  transformations  at 

random,  according  to  the  probabilities  p{  —  say  Tt  is  chosen,  and  apply  it  to  X0,  thereby  obtaining 

Xt  =  TtX0.  Then  choose  a  transformation  again  at  random,  independent  of  the  previous  choice,  and 

apply  it  to  Xlt  thereby  obtaining  X2.  Continue  in  this  fashion,  and  plot  the  orbit  {X„}.  The  result  is 

the  leaf  shown.  By  tabulating  the  frequencies  with  which  the  points  Xn  fall  into  the  various  pixels  of 

the  graphics  window,  one  can  actually  plot  the  empirical  distribution  — =-r  £  ,  using  a  grey  scale 

n+1 k=0  * 

to  convert  statistical  frequency  to  grey  level.  The  darker  portions  of  the  leaf  correspond  to  high 
probability  density. 


I.  Encoding  via  the  Collage  Property 

Given  a  target  set  C,  the  affine  mappings  T,  for  the  IFS  to  (nearly)  reproduce  it  are  obtained 
through  the  following  optimization  problem 

MAXIMIZE  [area^C  Cl  R^)  -  A  areaa(T,C\C)] 
affine  contractions  T,. 

where  R<  are  the  residual  sets  Rq  =  C, 


=  Ri-ATiC 
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and  A  is  an  adjustable  parameter.  This  optimization  problem  is  solved  sequentially  for  TltT2 _  The 

motivation  behind  this  technique,  and  an  adaptive  algorithm  developed  by  J.-P.  Vidal  for  solving  this 
problem  were  described  in  my  Second  Annua]  Report  (dated  12  April  1989),  and  in  my  Research 
Progress  and  Forecast  Report  for  the  First  Additional  Period  of  Research  (dated  14  October  1988). 

II.  Encoding  via  Convexity 

The  PI  has  shown  that  the  extreme  points  of  the  attractor  C  for  an  IFS  are  invariant  under  the 
inverse  dynamics.  That  is,  every  extreme  point  is  the  image  of  some  extreme  point  under  one  of  the 
afline  mappings.  Equivalently 


ECU  TjE  (1) 

l 

where  E  is  the  set  of  extreme  points  of  (the  convex  hull  of)  C,  and  T;  are  the  affine  maps  for  an  IFS 
which  generates  C.  This  is  illustrated  in  Figure  2.  Since  an  afline  map  T:  R*J  —►  Ra  is  completely 
determined  by  its  action  on  three  points,  this  condition  (1)  can  be  used  to  partially  or  fully  determine 
the  maps  Tt-  from  knowledge  of  the  image  C.  In  fact  condition  (1)  was  observed  by  the  PI  to  hold  in 
more  general  settings  of  Julia  sets,  where  the  maps  T,  are  the  inverse  branches  of  a  complex 
polynomial  f(z).  In  this  case  f(E)  C  E,  where  E  is  the  set  of  extreme  points  of  the  Julia  set. 

J.-P.  Vidal  implemented  an  algorithm  for  finding  all  possible  sets  of  affine  maps  T{  which  are 
consistent  with  (1).  This  was  described  in  my  First  Annual  Technical  Report  (dated  15  April  1988), 
and  in  my  Research  Progress  and  Forecast  Report  for  the  First  Additional  Period  of  Research  (dated 
14  October  1988).  It  was  published  in  [5].  (See  the  list  of  publications  at  the  end  of  this  report.) 

III.  Encoding  via  Discrete  Optimisation 

The  PI,  together  with  M.  Perrugia,  studied  a  discrete  IFS  model,  whereby  a  continuous  IFS 
process  is  generated,  and  every  point  along  the  orbit  is  sequentially  replaced  with  the  center  of  the  pixel 
in  which  it  lies.  It  was  shown  that  if  the  continuous  IFS  is  strictly  contractive,  then  as  pixel  sise  — *  0 
any  sequence  of  stationary  distributions  for  the  discrete  IFS  converges  weakly  to  the  (unique) 
stationary  distribution  for  the  continuous  IFS.  The  discrete  IFS  need  not  have  a  unique  stationary 
distribution,  and  conditions  ensuring  ergodicity  of  the  discrete  IFS  when  pixel  size  is  sufficiently  small 
are  the  subject  of  M.  Perrugia's  Ph.D.  thesis. 

The  discrete  IFS  leads  to  a  discrete  formulation  for  the  encoding  problem.  Given  the 
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stationary  distribution  x  satisfying  xP  =  x,  where  P  is  a  transition  probability  matrix  having 
precisely  N  non-zero  entries  in  each  row,  recover  P.  Typically  the  dimensions  of  P  (say 
(1024)3x(1024)3  for  a  high  resolution  screen)  are  much  larger  than  N  (which  is  on  the  order  of  10). 
There  is  in  fact  more  structure  on  account  of  the  afTme-ness,  so  that  the  non-zero  entries  of  P  are 
constrained  to  form  generalized  diagonals.  This  was  described  in  my  Second  Annual  Technical  Report 
(dated  12  April  1989). 

IV.  Encoding  via  the  Radon  Transform 

Given  a  random  variable  X  €  R2,  the  Radon  transform  R(A,z)  is  defined  by 

R(A,z)  =  P(X  ■  A  <  z), 

where  A  €  R2  is  a  unit  vector  and  z  €  R.  If  X  has  the  stationary  distribution  for  an  IFS  with  affine 
maps  Tjixi-*  A,x  +  b<  and  probability  weights  pi(  then  the  invariance  condition  can  be  expressed  as 

R(A,z)  =  ?  p,R(A,.,z,.)  (2) 

l 


where 


\  —  AjA  _  z  —  b-  A 

'  “  IIAJAH*  ~  HAjAH  ' 


(If  AjA  =  0  then  define 


R(V,)  =  *{b..A<z}  •) 


The  PI  has  been  developing  an  algorithm  for  recovering  the  maps  T,  and  weights  p<  via  (2), 
through  knowledge  of  R(A,z).  It  proceeds  by  finding  eigenvectors  of  the  matrix  parts  Ait  through 
invariant  half-planes.  In  fact  this  work  was  what  led  to  discovery  of  the  extreme  point  invariance 
condition  described  above  in  §11.  Indeed  one  can  use  (2)  to  analyze  supporting  hyperplanes  X  •  X1  <  z' 
of  the  image,  for  which  (A'.z')  will  be  on  the  boundary  of  the  support  of  R(A,z). 

There  is  a  similar  difference  equation  for  the  Fourier  transform  (i.e.,  characteristic  function) 

^(u)  =  EeiX  u,  u  €  R3. 
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Namely, 

y>(u)  =  ?  p ft**  Uy>(A,‘u).  (3) 

The  PI,  together  with  R.  Shenk  (GA  Tech.),  has  also  been  developing  an  algorithm  for  solving  the 
encoding  problem  based  on  (3). 

V.  Mixing  Images 

Mixing  theory  involves  an  MxM  transition  probability  matrix  P  =  (p4  •)  and  M2  (Borel) 
probability  measures  €  ^(G),  where  G  is  the  semi-group  of  affine  transformations  R2  — ►  R2.  One 
imagines  that  M  computer  screens  are  each  generating  a  stochastic  process  (Xn(i):  n  >  0),  for  1  <  i  < 
M.  (That  is,  i  indicates  screen  number.)  Let  (In(i):  n  >  1)  be  an  i.i.d.  sequence  of  “switching 
variables”  distributed  like  the  rows  of  P;  that  is, 

P(I»(i)=j)  =  P0. 


The  mixed  process  then  evolves  as 

X„+i(i)  =  g„+i(ij)Xn(j),  where  j  =  In+1(i). 

Here  (gn(ij):  n  >  1)  is  an  i.i.d.  sequence  distributed  like  Thus  if  one  represents  the  screens  as 
nodes  of  a  directed  graph,  with  processes  evolving  in  them,  then  the  (n  +  1)1'  point  generated  in  node  i 
i»  g»+i(ij)  applied  to  the  n<h  point  from  node  j.  The  index  In+1(i)  represents  the  crossover  from  node  j 
to  node  i,  along  edge  (j,i).  Figure  3  illustrates  a  typical  mixing  scenario.  This  structure  generates  the 
coupling  necessary  to  mix  the  individual  images  that  would  have  been  generated  in  each  screen  by  the 
maps  (gn(i,i):  n  >  1)  alone.  The  details  of  this  work  appear  in  the  Second  Annual  Technical  Report 
(dated  12  April  1989),  the  Research  Progress  and  Forecast  Report  for  the  First  Additional  eriod  of 
Research  (dated  14  October  1988),  the  First  Annual  Technical  Report  (dated  15  April  1988)  and  the 
Research  Progress  and  Forecast  Report  (dated  11  Sept.  1987).  It  was  published  in  [1],  [3],  [5]. 

VI:  Animation  and  Veetorisation 

Suppose  one  has  two  IFS,  each  corresponding  to  a  specific  image.  One  can  set  up  a  continuous 
flow  of  images  from  one  to  the  other  by  interpolating  the  IFS  parameters.  What  makes  IFS  animation 
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even  more  special  is  the  ease  with  which  one  can  rotate,  scale,  change  perspective  or  vantage  point, 
room  in  and  out,  or  perform  any  affine  transformation  on  an  image.  This  is  because  the 
transformations  can  be  built  directly  into  the  IFS. 

IFS  animation  is  highly  parallel.  The  images  for  the  various  frames  can  all  be  generated 
simultaneously,  since  there  are  no  dependencies  among  the  different  frames.  Furthermore  the  same 
sequence  of  random  numbers  can  be  used  for  all  the  images.  The  PI  was  able  to  exploit  the  vector  and 
multi-processor  capabilities  of  the  CRAY  Y-MP/832  at  the  Pittsburgh  Supercomputing  Center  to 
generate  real-time  animation.  This  was  described  in  the  Second  Annual  Technical  Report  (dated  12 
April  1989),  and  a  demo  video  tape  was  sent  in  to  the  Program  Director.  It  was  published  in  [2]. 

VII:  Wavelets  and  Curve  Generation 

Many  algorithms  for  generating  computer  images  today  involve  a  recursive  tree-traversal. 
These  include  wavelets  and  solutions  to  dilation  equations,  sub-division  refinement  methods  for 
generating  B-spIines  and  Bezier  curves  and  line  averaging  methods  for  interpolants.  Using  ergodic 
theory  the  PI  is  able  to  produce  random  algorithms,  in  a  very  general  setting,  which  generate  the  same 
images  as  the  recursive  ones.  These  images  become  attractors  of  random  dynamical  systems,  and 
evolve  simply  as  the  trajectory  of  a  single  orbit.  The  random  algorithms  are  very  fast,  involving  only 
affine  arithmetic,  and  are  efficiently  and  highly  parallelizable. 

In  sub-division  and  line  averaging  methods  for  curve  generation  one  is  given  mxm  row 
stochastic  matrices  P(0),  ...,  P(N  — 1)  (i.e.,  all  row  sums  are  one  —  negative  entries  allowed)  and  a  set 
of  “control  points”  vlt  ...,  vm  6  Rrf.  These  points  form  the  vertices  of  a  “control  polytope”  C  C  Rd, 
which  one  can  associate  with  the  mxd  matrix 


(This  correspondence  is  well-defined  if  one  thinks  of  C  as  an  “ordered”  polytope.)  If  P  =  (p4j)  is  an 
mxm  row  stochastic  matrix  then  one  can  identify  an  action  of  P  on  C;  namely,  PC.  Equivalently  PC 
is  the  polytope  with  vertices  v{,  ...,  vj„  where 

v!  =  ?p -X  • 

Under  suitable  conditions  on  the  P(i)’s  it  can  be  shown  that 
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P(in)  P(h)C 

converges  to  a  singleton  as  n  — ►  oo,  for  any  sequence  (i^ij,  ...).  By  choosing  the  P(i)’s  appropriately 
one  generates  a  smooth  curve  jC  (the  “attractor”)  as 

•A  =  U(lim  P(in)  •••  P(ii)C:  all  sequences  (ij,i2,  ...)) 

Let  Hm  C  Rm  be  the  hyper-plane 


Hm  —  {x  €  Rm:  ?  Xj  —  1}, 

i 

and  let  II  denote  the  projection  Rm  — .  Rm_1  which  lops  off  the  last  component 


■ 

- 

Xl 

*1 

H: 

1— ► 

... 

Xm 

«-H 

1 

£ 

X 

Given  any  mxd  matrix  C  there  is  a  unique  affine  transformation  T:  Rm-1  — ►  Rd  making  the  following 
diagram  commute. 


n 


Hm  - >  Rm_1 


Denote  T  =  T(C).  It  can  be  written  out  in  coordinate  form 

Tx  =  [Vj-Vm  |  •••  |  Vm_j-Vm]x  +  Vm 

where  v‘,  v„  are  the  row  vectors  of  C.  Similarly  if  P  is  an  mxm  row  stochastic  matrix  then  there 
is  a  unique  affine  transformation  T:  Rm_1  -*  Rm_1  making  the  following  diagram  commute. 
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In  fact  T  =  T(PH‘),  so  that  the  matrix  part  of  T  is  given  by  A  =  (pji  —  pmi)mT11  and  the 
translational  part  of  T  is  the  column  vector  (Pm*)™."1. 

The  Pi’s  random  algorithm  for  generating  the  curve  J.  follows.  The  inputs  are  the  row 
stochastic  matrices  P(i)  and  the  control  polygon  C.  Define  T(i)  from  P(i)  via  the  commutative 
diagram  above. 

Algorithm  (Sub-Division  Methods): 
initialize  x  =  (the  fixed  point  of  T(0))  €  Rm_1 
for  n  =  1,  10000 
plot  T(C)x  €  H* 

choose  i  €  {0,  ...,  N  — 1}  randomly 
x  -  T(i)x 

This  algorithm  is  illustrated  in  Figure  4.  The  curves  in  Figs.  4a)-e)  are  all  quadratic  B-splines. 
They  correspond  to  different  placements  of  the  control  points  v{.  Observe  how  this  algorithm 
parallelizes.  Instead  of  running  an  orbit  of  length  10000  on  one  processor,  one  can  run  (say)  10 
processors,  each  generating  only  1000  points  of  an  orbit. 

Wavelets  are  compactly  supported  functions  W(x),  x  €  R,  having  the  special  feature  that  the 
family  {W(2;x— k):  j  >  0,  k  €  Z)  of  translates  and  dilations  are  orthogonal  (in  the  Lj-sense).  These 
functions  W(2;x— k)  are  local  in  both  space  and  frequency,  and  thus  afford  a  new  type  of  basis  for 
orthogonal  expansions. 

Daubechies’  wavelets  are  constructed  through  solutions  g(x),  x  €  R,  of  certain  dilation 
equations 

rg(x)  =  Ea*g(2x-k) 

{  Eg(x— k)  =  l 


(4) 
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Precisely,  W  is  given  in  terms  of  g  as 


W(x)  =  £(-l)V*g(2*— k)  (5) 

If  the  coefficients  at  are  chosen  in  a  very  special  way,  then  W  will  have  the  desired  orthogonality 
property.  Precisely  at  =  0  for  k  <  0  or  k  >  m,  where  m  =  2p —  1;  and  the  2p  coefficients  a0,  ...,  am 
have  to  satisfy  the  2p  conditions 


£*2*  —  1*  £*2*+i  —  1 

£(-l)Vat  =  0;  t  =  1 . P-1 

Ea*a**2«  =  °;  1  =  1 . P"1 

Then  W  =  W p  given  by  (5)  will  have  the  desired  orthogonality  property. 

Wavelet  generation  can  be  put  into  the  preceding  framework  by  taking  N  =  2,  and  setting 


P(0)  =  (a  P(D  =  (<h,-i)Z=v 


and  taking  as  control  points 


(6) 


=  tj  e  Rm-\  1  <  j  <  m— 1, 
where  e;-  are  the  standard  unit  vectors,  and 


vm  =  0. 

Define  T(0),  T(l)  in  terms  of  P(0),  P(l)  via  the  commutative  diagram  above.  Then  the  Pi’s  random 
algorithm  for  generating  the  solution  g  of  the  dilation  equations  is  as  follows. 

Algorithm  (Solution  of  Dilation  Equations): 
initialise  a  =  0,  x  =  (the  fixed  point  of  T(0))  €  Rm_1 
for  n  =  1,  10000 


k=l 


plot  (a,x,),  ...,  (a+m-1,  xm) 
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choose  i  €  {0,1}  randomly 
a  —  j(a+i) 
x  «—  T(i)x 

This  is  illustrated  in  Figures  5  and  6.  The  corresponding  algorithm  for  generating  wavelets, 
based  on  (5),  is  as  follows. 

Algorithm  (Daubechies’  Wavelets): 

initialize  a  =  0,  x  =  (the  fixed  point  of  T(0))  €  Rm_1 

for  n  =  1,  10000 

m— 1 

Xrn  =  1  -  J2  xt 
k=l 

min(m,m— £) 

plot  (a  +  £,  £  (— l)*-  a.**  £ ),  — m  +  1  <  l  <  m 

V  1  k=max(0,l-£)  ' 

choose  i  €  {0,1}  randomly 

*  —  J(®  +  5) 

x  -  T(i)x 

This  is  illustrated  in  Figures  7  and  8. 

In  all  of  these  algorithms  it  is  not  necessary  to  initialize  x  as  the  fixed  point  of  T(0).  Rather 
one  can  initialize  x  arbitrarily  (say  x  =  0  €  Rm_I),  and  skip  the  plotting  until  n  =  10  or  20.  The 
point  is  that  no  matter  how  x  is  initialized,  the  IFS  orbit  will  quickly  move  into  the  attractor.  This 
work  was  published  in  [6],  [7],  [8]. 

What  Next? 

There  are  many  important  follow-up  problems  to  attack  and  algorithms  to  develop  for  IFS 
theory.  Listed  below  are  some  objectives  the  PI  has  already  begun  working  on. 

1.  Conditions  ensuring  ergodicity:  When  generating  curves  using  sub-division  methods  one  is 
given  mxm  row  stochastic  matrices  P(0),  ...,  P(N— 1)  as  described  above.  If  these  matrices  have  all 
non-negative  entries  then  under  an  irreducibility-type  condition  it  can  be  shown  that  the  random 
algorithm  above  for  sub-division  methods  converges.  More  precisely,  it  can  be  shown  that  the 
underlying  Markov  chain  is  ergodic.  In  fact  it  suffices  to  show  that 


fi™  h  logllA(ii)  •••  A(i„)||  <  0 


(7) 


for  any  sequence  —)i  where  A(i)  is  the  (m— l)x(m— 1)  matrix  part  of  the  affine  map  T(i) 

constructed  above  (via  the  commutative  diagram). 

On  the  other  hand  if  the  matrices  P(i)  are  allowed  to  have  some  negative  entries,  then 
conditions  ensuring  (7)  are  not  as  apparent.  In  particular  for  line  averaging  schemes  where  P(0),P(1) 
are  given  by  (6),  what  conditions  on  the  coefficients  at  ensure  (7)?  Or  if  the  coefficients  at  depend  on 
some  parameters  6 ,  what  are  the  critical  values  of  0  for  which  the  underlying  Markov  chain  switches 
from  recurrence  to  transience?  This  is  the  most  significant  and  pressing  question.  Partial  answers  are 
given  in  [6,  Thm.  II],  where  conditions  are  presented  under  which  there  will  exist  some  operator  norm 
making  A(0),  ...,  A(N  —  1)  strict  contractions.  But  more  work  needs  to  be  done  on  this  question. 

2.  Surface  generation:  The  work  on  curve  generation  for  sub-division  methods  and  wavelets 
needs  to  be  extended  to  surfaces.  Surface  generation,  like  curve  generation,  can  be  set  up  through 
refinement  schemes;  and  the  Pi’s  technique  produces  random  algorithms  which  generate  the  same 
surfaces  as  the  recursive  algorithms.  This  needs  to  be  carried  out. 

3.  Dimension  analysis:  The  fractal  dimensions  of  the  curves  and  surfaces  produced  by  the 
above  random  algorithms  need  to  be  determined.  Fractal  dimension  is  a  major  characteristic  of  surface 
roughness.  In  special  cases  the  dimension  of  an  IFS  attractor  is  known,  and  this  analysis  needs  to  be 
extended  to  the  setting  of  sub-division  methods  described  above,  involving  row  stochastic  matrices 
P(0),  ...,  P(N— 1).  Knowing  the  dimension  can  also  lead  to  determination  of  critical  parameter  values, 
as  described  in  1.  above. 

4.  Extreme  points  of  attractors:  Let  J  be  the  Julia  set  for  a  complex  polynomial  f(z),  and  let 
E  be  the  set  of  extreme  points  of  (the  convex  hull  of)  J.  As  mentioned  in  §11  above,  the  PI  has 
observed  in  some  computer-graphical  experiments  that  f(E)  C  E.  This  needs  to  be  investigated,  to  see 
whether  this  invariance  property  applies  to  more  general  dynamical  systems  than  affine  IFS. 

5.  Animation  encoding:  As  a  natural  follow-up  to  the  work  on  fractal  animation  described 
above  in  §VI,  the  PI  has  been  looking  into  encoding  for  animation.  Here  the  data  compression  ratios 
are  enormous,  since  the  encoding  of  two  “endpoint*  images  suffices  to  generate  the  intermediates. 
Furthermore  in  certain  respects  animation  encoding  is  easier  than  still  image  encoding.  This  is  because 
a  dynamic  sequence  of  images  often  conveys  additional  information  about  the  individual  still  images, 
such  as  segmentation.  Velocity  tracking  of  boundaries  and  key  features  of  an  image  can  be  used  to 
decide  where  to  position  the  temporal  IFS  linear  interpolation  points  (i.e.,  where  to  break  the 


animation  up  into  “piecewise  linear”  video  segments),  and  how  to  identify  the  images  as  1FS  mixtures. 
Animation  is  potentially  the  most  exciting  application  of  IFS  encoding. 
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Tx(x) 

T3(x) 


°o  am  a)  «•> 

0.355  -0.355  \  ( 0.266  ^  T(  . 
0.355  0.355  j  *  +  \  0.078  )  4{X) 


0.5  0  \  /  0.25  \ 

0  0.5  /  +  \  0.4  J 

0.355  0.355  ^  /  0.378  \ 

-0.355  0.355  J  *  +  \  0.434  ) 


P\  =  0.5  p2  =  0.168  P3  s  0.166  Pi  =  0.166 


Figure  1:  Maple  Leaf 

This  image  waa  generated  using  the  basic  IFS  algorithm.  The  affine  maps  listed  above  were  composed 
randomly,  according  to  the  respective  probabilities  shown  beneath  them.  The  darker  pixels  had  a 
higher  proportion  of  points  in  the  orbit  fall  into  them.  (The  window  here  is  0  <  x,y  <  1.) 
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/  0.5  0.5  \  (  0.125  \  T(.  (  0.5  0.5  \  /  -0.125  \ 

V  -0-S  0.5  JI+(  0.625  )  ™  =  {  -0.5  0.5  )  *  +  (  0.375  J 


Pi  =  0.5  pi  —  0.5 


jt  mm*  jt 


Figure  2:  Convex  Hull  of  an  Attractor 

Each  extreme  point  is  the  image  of  some  extreme  point  under  T,  or  T2.  The  interior  angles  of  this 
convex  hull  are  all  135*.  This  is  because  the  maps  T,  and  Ta  are  conformal.  (The  window  here  is 
— 1  <  x,y  <  1.) 
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read  vt,  v3,  v3,  v4 
initialize  x  =  y  =  z  =  0 
for  n  =  1,  10000 

plot  K-v4  I  v3-v4  I  v3-v4]x  +  v4 
choose  i  6  {0,1}  randomly 
if  i  =  0 

then  x  «—  x  +  ^y 

J- + Jy  + 

z  •—  1  —  x  —  y 
else  z^-jx  +  ly  +  gz 

y - h*+ 4y 

x  < —  0 


Figure  4:  Random  Algorithm  to  Generate  Quadratic  B-Splines 
Pick  any  points  v1(  v2,  v3,  v4  and  this  random  algorithm  produces  their  quadratic  spline  fit.  This 
curve  passes  through  vlt  ^  v2  +  ^  v3,  v4  and  is  tangent  to  the  lines  VjVj,  v2v3,  v3v4  respectively,  at 
these  points.  The  row  stochastic  matrices  here  are  given  by 


1 

0 

0 

0 

0 

1/2 

1/2 

0 

1/2 

1/2 

0 

0 

.  P(l)  = 

0 

1/4 

3/4 

0 

0 

3/4 

1/4 

0 

0 

0 

1/2 

1/2 

0 

1/2 

1/2 

0 

0 

0 

0 

1 

This  algorithm  can  be  optimized  by  running  it  in  parallel,  breaking  up  the  loop  above  into  many 
smaller  loops. 


18 


initialize  a  =  x  =  0;  y  =  0.6154 
for  n  =  1,  10000 
z  =  1  —  x  —  y 
plot  (a*x),  (a+1,  y),  (a+2,  z) 
choose  i  €  {0,1}  randomly 
if  i  =  0 
then  a  «—  0.5a 

y  «—  0.8  —  0.6x  -  0.3y 
x  «—  0.8x 

else  a  «—  0.5a  +  0.5 
x  <—  0.5x  •+■  0.8y 
y  *-  0.5  -  0.3y 

Figure  5:  Random  Algorithm  to  Generate  the  Solution  of  the  Dilation  Equations 
This  random  algorithm  plots  the  solution  of  the  dilation  equations  (4)  with  (non-zero)  coefficients 

a0  =  ®i  =  0.5,  s  0.2,  83  ~  0.5 


1 


2 


3 
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initialize  a  =  x  =  0;  y  =  1.366 
for  n  =  1,  10000 
z  =  1  —  x  —  y 
plot  (a,x),  (a+1,  y),  (a+2,  z) 
choose  i  €  {0,1}  randomly 
if  i  =  0 
then  a  «—  0.5a 

y  - - 0.366x  +  0.5y  +  0.683 

x  -  0.683x 
else  a  «—  0.5a  +  0.5 
z  «—  x 

x  «-  1.183x  +  0.683y 
y  < - 1.366z  -  0.866y  +1.183 


Figure  6:  Random  Algorithm  to  Generate  the  Solution  of  the  Dilation  Equations 
This  random  algorithm  plots  the  solution  of  the  dilation  equations  (4),  used  to  construct  Daubechies’ 
wavelet  W4.  The  (non-zero)  coefficients  are 


»  _  1  +  43 

—  4  i 


.  -3  +  43 

ai  —  5  i 


a  _3-43 

—  4  i 


a  -L=JI 

a3  —  4 
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initialize  a  =  x  =  0;  y  =  1.366 
for  n  =  1,  10000 

bt  =  — 0.183x;  b2  =  -0.5x  -  0.183y;  b3  =  0.866x  -  0.317y  -  0.183 
b4  =  0.5x  +  1.183y  -  0.5;  b5  =  0.683(l-x-y) 

plot  (a— 1,  bj,  (a— 0.5,  b2— bt),  (a,  b3-b2),  (a+0.5,  b4  — b3),  (a+1,  b5  — b4),  (a+1.5,  — b5) 
choose  i  €  {0,1}  randomly 
if  i  =  0 
then  a  «—  0.5a 

y  - - 0.366x  +  0.5y  +  0.683 

x  v-  0.683x 
else  a  ♦—  0.5a  +  0.25 
z  <—  x 

x  <-  1.183x  +  0.683y 
y  - - 1.366z  -  0.866y  +  1.183 

Figure  7:  Random  Algorithm  to  Construct  Daubechiea’  Wavelet  W4 
This  wavelet  is  constructed  out  of  the  curve  from  Fig.  6  above,  via  (5).  There  is  no  need  for  a 
recursive  tree  algorithm  —  ergodic  theory  does  all  the  running  around. 
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initialize  r  =  z0  =  z8  =  xt  =  0;  x2  =  1.286;  x3  =  —0.386;  x4  =  0.095 

»o  =  0.470;  at  =  1.141;  a2  =  0.650;  a3  =  -0.191;  a4  =  -0.121;  a5  =  0.050 
for  n  =  1,  10000 

X8  =  1  —  xx  -  x2  -  x3  -  x4;  a„  < - a0;  a2  < - a2;  a4  < - a4 

for  l  =  —4,  5 

mj  =  max(0,  1  —  i),  m2  =  min(5,  5  —  6);  y  =  0 
Fork  =  m1,m2 

y  -  y  + 

endfor 

plot  (r  +  | ,  y) 
endfor 

choose  i  €  {0,1}  randomly 

r  «—  0.5r  +  0.25i;  a„  < - a0;  a2  - a2;  a4  < - a4 

for  t  =  1,  5 

*<  -  xv  xi  *-  0 

endfor 

for  €  =  1,4 

mx  =  max(0,  2€ — 6);  m2  =  min(5,  2€  —  1);  j  =  2<  +  i  —  1 
for  k  =  mu  m2 

xt*~xt  +  *kzj-k 
endfor 

endfor 

endfor 


Figure  8:  Random  Algorithm  to  Construct  Daubechies’  Wavelet  Wa 
The  (non-zero)  coefficients  at  are  given  by 
ag  =  0.470,  at  =  1.141,  a2  =  0.650,  a3  =  —0.191,  a4  =  —0.121,  a5  =  0.050 


Figure  8 
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Published  Articles  under  this  Grant 

1.  Barnsley,  M.  F.,  Berger,  M.  A.  and  Soner,  H.  M.,  Mixing  Markov  chains  and  their  images,  Prob. 
Eng.  Inf.  Set.  2  (1988),  387-414. 

2.  Barnsley,  M.  F.  and  Berger,  M.  A.,  Pictures  worth  a  million  words,  in  Projects  in  Scientific 
Computing:  1988-1989,  Pittsburgh  Supercomputing  Center,  34-35. 

3.  Berger,  M.  A.,  Encoding  images  through  transition  probabilities,  Math!.  Compvt.  Modelling  11 
(1988),  575-577. 

4.  Berger,  M.  A.,  Images  generated  by  orbits  of  2-D  Markov  chains,  CHANCE  2  (1989),  18-28. 

5.  Berger,  M.  A.,  Random  affine  IFS:  mixing  and  encoding,  submitted  (Trans.  Amer.  Math.  Soc.). 

6.  Berger,  M.  A.,  Random  affine  IFS:  smooth  curve  generation,  submitted  (SIAM  Review). 

7.  Berger,  M.  A.,  Ergodic  theory,  tree  traversal  and  computer  image  generation,  submitted  (Proc. 
Nat’l  Acad.  Sci.) 

8.  Berger,  M.  A.,  Wavelets  as  attractors  of  uynamical  systems,  submitted  (Bull.  Amer.  Math.  Soc.) 

9.  Berger,  M.  A.,  Review  of  Fractals  Everywhere,  Stochastics  and  Stochastic  Reports,  in  press. 

10.  Berger,  M.  A.  and  Soner,  H.  M.,  Random  walks  generated  by  affine  mappings,  J.  Theor.  Prob.  1 

(1988),  239-254. 

Article  [2]  above  was  prepared  by  science  editor  Michael  Schneider  (Pittsburgh  Supercomputing 
Center). 
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Interactions  under  this  Grant 

The  PI  gave  the  following  presentations: 

1.  Sixth  International  Conf.  on  Math.  Modelling,  held  at  Washington  Univ.,  Aug.  10-14,  1987  (host: 
Dr.  I.  Rodin). 

2.  Invited  talk  for  the  seminar  run  by  the  Pittsburgh  Supercomputing  Center  in  October,  1987  (host: 
Dr.  R.  Roskies). 

3.  AFOSR  at  Bolling  Air  Force  Base  on  Feb.  24,  1988  (host:  Dr.  A.  Nachman). 

4.  Invited  talk  in  Michael  Barnsley’s  minisymposium  on  chaos  at  the  annual  SIAM  meeting  in 
Minneapolis,  July  10-15,  1988  (host:  Dr.  M.  Barnsley). 

5.  Invited  talk  for  the  seminar  run  by  the  image  processing  group  (Grenander,  McClure,  Geman  and 
Gidas)  in  the  Division  of  Applied  Mathematics  at  Brown  University,  Feb.  15,  1989  (host:  Dr.  B. 
Gidas). 

6.  NIST  in  Gathersburg,  MD  on  March  29,  1989  (host:  Dr.  F.  Sullivan). 

7.  NSF  in  Washington,  DC  on  March  30,  1989  (host:  Dr.  R.  Chin). 

8.  Special  3-day  Lecture  Series  at  Allegheny  College,  April  10-12,  1989. 

9.  Invited  Address,  SIAM  Conference  on  Dynamical  Systems,  May  7-10,  1990  (host:  Dr.  H.  Stech). 

10.  Invited  talk  at  the  meeting  on  Lyapunov  Exponents,  Mathematisches  Forschungsinstitut 
Oberwolfach  (Black  Forest,  Federal  Republic  of  Germany),  May  27-June  2,  1990  (host:  Dr.  L. 
Arnold,  Bremen). 

11.  Invited  talk  at  the  First  IFIP  Conference  on  Fractals,  FRACTAL  90,  Lisbon,  Portugal,  June  6-8, 
1990  (host:  Dr.  H.-O.  Peitgen,  Bremen). 
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Participating  Professionals 

Some  of  the  Pi’s  research  on  image  encoding  was  carried  out  with 

(1)  Michael  Barnsley  (Prof,  of  Math.,  GA  Tech.); 

(2)  Bill  Eddy  (Prof,  of  Statistics,  CMU); 

(3)  John  Elton  (Prof,  of  Math.,  GA  Tech.); 

(4)  Jeff  Geronimo  (Prof,  of  Math.,  GA  Tech.); 

(5)  Mario  Perrugia  (Ph.D.  student,  Statistics,  CMU  —  partially  funded  by  AFOSR); 

(6)  Ron  Shonkwiler  (Prof,  of  Math.,  GA  Tech.); 
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